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We examine the effective counterion-mediated electrostatic interaction between two like-charged 
dielectric cylinders immersed in a continuous dielectric medium containing neutralizing mobile coun- 
terions. We focus on the effects of image charges induced as a result of the dielectric mismatch 
between the cylindrical cores and the surrounding dielectric medium and investigate the counterion- 
mediated electrostatic interaction between the cylinders in both limits of weak and strong elec- 
trostatic couplings (corresponding, e.g., to systems with monovalent and multivalent counterions, 
respectively). The results are compared with extensive Monte-Carlo simulations exhibiting good 
agreement with the limiting weak and strong coupling results in their respective regime of validity. 
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I. INTRODUCTION 

Electrostatic interactions of charged macromoleculcs 
and colloids are often governed by small oppositely 
charged ions (counterions) that maintain global elec- 
troneutrality. These counterion-mediated interactions 
play a fundamental role in classical charged (Coulomb) 
fluids that are abundant in biological and soft matter 
context [l|, 0] and include many charged macromolecules 
(such as nucleic acids DNA and RNA, actin filaments, 
microtubules and globular proteins) , affecting their func- 
tional, structural and dynamical behavior. In spite of 
the importance of electrostatic interactions, there is no 
general method that would allow for an accurate pre- 
diction of electrostatic effects in all regions of the pa- 
rameter space, defined by the surface charge density 
of macroions, charge valency of counterions, dielectric 
mismatches between the often hydrophobic core of the 
macromolecule and the surrounding aqueous medium, 
etc. Often the electrostatic interactions are treated on 
the Poisson-Boltzmann (PB) level leading to effective in- 
teractions which turn out to be always repulsive between 
like-charged macromoleculcs. Conceptually, the PB ap- 
proach corresponds to a mean-field treatment of electro- 
static interactions and is asymptotically valid only for 
sufficiently large separations between macromolecules, 
low enough surface charge densities and low counterion 
valency It characterizes a situation where electro- 
static fluctuations and correlations due to the counteri- 
ons are negligible. There are other regions in the pa- 
rameter space of charged macromoleculcs where one ex- 
pects the mean-field framework to break down leading 
to the emergence of a completely different non-PB-type 
physics. A notorious example is the phenomenon of like- 
charge attraction, which emerges between highly charged 
macroions or in the presence of high valency counterions 
and has been at the focus of both experimental [3l-fl0j 
and theoretical investigations in recent years (see Refs. 
0, [Tl| - |2fjj for an extensive reference list) . 

It appears to us that among the most important re- 



cent advances in this field has been the systematization 
of non-PB effects based on the notions of weak coupling 
(WC) and strong coupling (SC) approximations. These 
terms refer to the strength of electrostatic coupling in the 
system and may be understood conceptually in terms of 
the electrostatic interactions of mobile counterions with 
fixed external charges (macroions) in the system when 
compared with direct electrostatic interactions between 
the counterions themselves. This latter contribution may 
be characterized in terms of the Bjerrum length, 



e 2 /(4iree k B T), 



(1) 



which corresponds to the separation at which two unit 
charges, eo, interact with thermal energy k B T in a 
medium of dielectric constant s (in water and at room 
temperature, the value is £q « 0.7 nm). If the charge of 
the counterions is +qeo then the Bjerrum length scales 
as q 2 £B- The interaction of counterions with macroion 
charges (of surface charge density —a s ) can be charac- 
terized in terms of the so-called Gouy- Chapman length, 



(i = e /(2Trq£ B o- s ), 



(2) 



which gives the separation at which the counterion- 
surface interaction energy equals k B T. The ratio of these 
two fundamental length scales introduces a dimensionlcss 
parameter 



a = q 
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(3) 



which is known as the (Netz-Moreira) electrostatic cou- 
pling parameter S [3 and quantifies the strength of 
electrostatic coupling in the system. This parameter is 
closely related to the plasma parameter of ionic systems 
[2ll ] and may be written also in terms of the typical lat- 
eral spacing, a±, between counterions in the proximity 
of a charged surface, i.e., a±/ fi ~ vS as set by the local 
electroneutrality condition a\ ~ qeo/a s . 

It then follows that in the weak coupling regime, de- 
fined by S <§; 1, the width of the Gouy- Chapman layer 
11 is much larger than the separation between two neigh- 
bouring counterions in the counterion layer and thus it 
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behaves basically as a three-dimensional gas. Each coun- 
terion in this case interacts with many others and the 
collective mean-field approach of the Poisson-Boltzmann 
type is thoroughly justified. On the other hand, in the 
strong coupling regime, defined by S 3> 1, the mean dis- 
tance between counterions, a±, exceeds the layer width 
and thus the counterion layer behaves essentially as a 
two-dimensional layer [ljj]. In this case the mean-field 
approach breaks down as each counterion is isolated lat- 
erally in a relatively large correlation hole of size aj_. 
In fact, as each counterion can move almost indepen- 
dently from the others along the direction perpendicular 
to the charged surface, the properties of the system are 
dominated by single-particle contributions on the lead- 
ing order, which is in stark contrast with the collective 
mean-field picture and emerges as a direct consequence 
of strong electrostatic correlations in the system. The 
two dychotomous limits are characterised by a low (high) 
valency of the counterions, a small (large) value of the 
surface charge density and/or large (small) medium di- 
electric constant. 

Conceptually the study of the SC regime has been pio- 
neered in several recent works jFp - ll8l ] using various ana- 
lytical methods. It was shown [l4[ that both the WC and 
the SC limits may be described analytically as two ex- 
act limiting laws from a single general theory for classical 
Coulomb fluids: while the PB theory follows in the limit 
of 3 — > 0, a limiting single-particle SC theory follows 
in the limit S — t 00, which thus allow for an exact sta- 
tistical mechanical treatment of charged systems at two 
opposed limiting conditions. The parameter space in be- 
tween can be analyzed by approximate methods [22M24| , 
beingaccessible effectively only via computer simulations 
til [l7|, [H, Hi-iU. Exact solutions for the whole 
range of coupling parameters are available only in one 
dimension 13211 . The WC-SC par adig m has been tested 
extensively [13, H El H-H M MM, 

and was found 

to describe computer simulations quantitatively correctly 
in the respective regimes of validity, thus providing a uni- 
fying view of the behavior of Coulomb fluids. The main 
facets of these results are retained even when the model 
system is generalized in order to include more realistic 
features describing the bathing solution or the nature of 
the fixed or mobile charges in the system. Though, for 
instance, multipolar charge distribution of mobile coun- 
terions 13311 or statistically disordered distribution of fixed 
charges [34| unavoidably introduce novel features in the 
counterion-mediated electrostatic interaction, the appli- 
cation of the same general philosophy embodied in the 
weak and strong coupling limits remains sound and valid. 
The only case where it needs to be amended in an es- 
sential manner is when the bathing solution contains a 
mixture of univalent as well as polyvalent salts, which 
incidentally are also the most common experimental con- 
ditions. In that case a more sophisticated mixed weak- 
strong coupling analysis is in order leading to qualita- 
tively different results [H HH . 

Though originally formulated in the context of pla- 



nar macroions, these advances have transpired also in 
the DNA-like models of polyelectrolytes which deal with 
electrostatic interactions in the context of macroions with 
cylindrical or indeed helical fixed charge distribution. In- 
deed this particular variant of counterion-mediated inter- 
actions has a venerable history [38l - l42l |. During the last 
two decades several distinct analytical approaches aimed 
in different directions improved the classical PB results 
for simple DNA-like models and revealed the importance 
of correlation effects as well as several other factors in- 
cluding the discrete or helical charge distributions, chain 
flexibility, finite molecular size, and dielectric inhomo- 
geneities [30|, HH, IUHfiSl • Guided by these developments 
we set ourselves the goal of systematically analyzing the 
interactions between cylindrical macroions mediated by 
mobile counterions in the presence of dielectric inhomo- 
gencities. 

Contrary to the case of planar macroions which can be 
completely characterized by a single dimensionless cou- 
pling parameter, cylindrical macroions require in general 
two dimensionless coupling parameters that consistently 
describe the range of validity of the strong and weak cou- 
pling approximations. The existence of two coupling pa- 
rameters is due to the simple fact that, if compared to the 
planar case, the cylinder has a finite radius that needs to 
enter the fundamental description of the problem. This 
other dimensionless parameter brought fourth by cylin- 
drical geometry of the macroion is nothing but the so- 
called Manning parameter [69j . For a cylinder of radius 
a and linear charge density A, it is given by 

£ = ^b- = ~ (4) 
e /i 

The Manning parameter thus represents the dimension- 
less linear charge density or, on the other hand, also 
the rescaled radius of the cylindrical charge distribution. 
Note that the two parameters (i.e., S and £) in the cylin- 
drical geometry arc independent and can be set sepa- 
rately. For double-stranded DNA (A rj 6eo/nm), the 
Manning parameter and the (Netz-Moreira) electrostatic 
coupling parameter are given by £ ~ 4.1 q and S ~ 2.8 q 3 
(at room temperature in water) in terms of the counte- 
rion valency q. 

In this paper we will analyze systematically the in- 
teractions between two cylindrical macroions character- 
ized by a fixed surface charge density as well as a fi- 
nite dielectric jump between the external dielectric back- 
ground (corresponding to a continuum solvent) and the 
hydrophobic cores of cylindrical macroions. Analytically 
derived results in both limits of strong and weak coupling 
will be compared directly with Monte-Carlo simulation 
results. The organization of this paper is as follows: We 
first introduce our model and study the image-charge ef- 
fects within the problem of a single charged cylinder with 
neutralizing counterions in Section [TTJ We then focus on 
the interaction between two identical charged dielectric 
cylinders within the WC and the SC theory as well as MC 
simulations in Section Mil The results in the presence 
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and absence of dielectric discontinuity at the cylindrical 
boundaries are analyzed in detail in order to bring up 
the effects due to the image charges in the two-cylinder 
geometry. 



II. ONE CHARGED CYLINDER 

Let us focus first on the problem of a single infinitely 
long uniformly charged cylindrical macroion of radius a. 
The charge of the cylinder is assumed to be distributed 
uniformly on its surface according to the charge distribu- 
tion function 



a{r) = -^-a dip - a) > 



(5) 



with A being the absolute linear charge density and p 
the radial coordinate from the cylinder axis. The cylin- 
der is arbitrarily chosen to be negatively charged and its 
electrical charge is thus neutralized by positively charged 
mobile counterions of valency +q, which are present in 
the region p > a only. 
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FIG. 1: Schematic top view of a single charged dielectric cylin- 
der of radius a and dielectric constant e , along with mobile 
neutralizing counterions of charge valency +q dispersed in a 
continuum solvent of dielectric constant e. The system is con- 
fined coaxially in a cylindrical box of radius a out (not shown) . 

In general, the interior of the cylindrical macro- 
molecule can be characterized by a different dielectric 
constant e' than the surrounding continuum solvent 
medium, e. This is certainly the case for DNA that has 
a hydrophobic inner core composed mostly of stacked ni- 
trogen bases that have a vastly different dielectric re- 
sponse from an aqueous solution [H| . Other charged 
polymers that usually do not share structural features 
with DNA can nevertheless also posses hydrophobic in- 
ner cores with a local dielectric response that differs from 
the one of the solvent. The hydrophobic apolar core of 
the macromolecular backbone would then have a static 
dielectric constant e' ~ 2 (hydrocarbon), while the sur- 
rounding continuum dielectric medium which is usually 
water, a polar associated liquid, would be characterized 
by e ~ 80. One can thus introduce the dielectric discon- 
tinuity parameter 



which measures the relative dielectric mismatch at the 
interface of the two materials. In a dielectrically homoge- 
neous system, we have A = and no image charge effects 
are present, while, for instance, in water-hydrocarbon 
systems, one has A = 0.95, which suggests strong image 
charge effects (note that |A| < 1 and the largest value 
for A is obtained when one medium is ideally polarizable, 
i.e., is an ideal metal). Therefore, a single cylinder can 
be described by three different dimensionlcss parameters, 
namely, the electrostatic coupling parameter S, the Man- 
ning parameter £ and the dielectric discontinuity param- 
eter A as defined in Eqs. ([3]), dU and ([5]), respectively. 

The presence of a dielectric inhomogeneity across the 
boundary of the cylinder, see Fig. [IJ influences the elec- 
trostatic potential that can be thus described by the 
Green's function connecting two points r, r' outside the 
cylinder as 



u(r, r') = u (r, r') + u im (r, r'), 



(7) 



where uq is the direct standard Coulomb kernel in the 
absence of dielectric inhomogeneities, 



u (r, r') 



f 



47T££o r — r' 



(8) 



and Uim is the "image correction" due to the dielectric 
jump just as in the case of a planar discontinuity. Un- 
fortunately in cylindrical geometry the concept of Kelvin 
image charges, that plays such a fundamental role for pla- 
nar dielectric discontinuities, is a bit diluted since the im- 
age correction can not be in general formulated in a way 
that would entail a summation of properly positioned 
point image charges. We nevertheless consistently refer 
to the effects of dielectric inhomogeneities in this system 
as dielectric image effects. If explicitly stated we believe 
this inconsistency in nomenclature can not be a source 
of confusion. 

Instead of using the concept of point image charges 
one must in fact explicitly solve the Poisson equation 
in cylindrical geometry specified by cylindrical coordi- 
nates (p, tp, z) with appropriate boundary conditions at 
the surface of the cylinder, which requires the normal 
component of the electric field to fulfill the following re- 
lation 



where 



££oE n \p =a + — e £ oEn\o=a a s 



X/(2wa), 



(9) 



(10) 



is the cylinder surface charge density. Solving the Poisson 
equation with the proper eigenfunction expansion (70j 
one can express the image part of the Green's function 
as 



u im (r,r') = 



f 



(6) 



x cos(fcAz) cos(mAip), 



dk £ m (ka)K m (kp')K m (kp) 
(11) 
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with the following definition 



2(2 - S m0 )AI„ 



(x) 



l + A 



3j/ m _)_i (i)+m/ m (x) 



2A/f m (x)' 



(12) 



where Ay> and Az are the angle and height differences be- 
tween the position vectors r and r', respectively. I m {x) 
and K m (x) are the modified Bessel functions of the first 
and the second kind that enter the eigenfunction expan- 
sion [zO]. 

Having derived the appropriate Green's function for 
the cylindrical geometry, we now investigate the behavior 
of counterions in proximity to a charged dielectric cylin- 
der in the two limits specified by the weak and strong 
coupling approximations that can be treated analytically. 
The behavior at intermediate coupling strengths will be 
analyzed by using extensive MC simulations. 



A. Weak coupling limit: Mean-field theory 

As noted before, the WC regime is characterized by 
a small coupling parameter S < 1, which is adequate 
for low valency counterions, small surface charge, high 
temperature or high solvent dielectric constant. In the 
strict limit of S — > 0, ignoring the usually small fluctua- 
tions around the mean field configuration, the system is 
described exactly by the mean-field PB theory [lfjfl . For- 
mally the PB equation corresponds to the saddle-point 
condition imposed on the action field-functional in this 
limit [tTJ] . For the present system, the PB equation gov- 
erning the mean electrostatic potential takes the form 



em 



' C e-PwK*) p>a , 



(13) 







P < a, 



where the right-hand side is obviously the number density 
of the counterions, i.e. 



i(r) =Ce^ eo9 ^ r) , p>a. 



(14) 



Note that due to axial symmetry and translational in- 
variance in the z direction, the problem is reduced to 
a one-dimensional formulation with V 2 = V p 2 and the 
solution of the PB equation depends only on the radial 
coordinate p. The analytical solution of the PB equa- 
tion in this case is well known and was first discussed in 
Refs. [39l - l4ll | within the so-called cell model where the 
cylindrical macroion is confined (coaxially) in an outer 
cylindrical boundary of radius a ou t ■ The cell model guar- 
antees a finite solution for the counterion density despite 
the extremely slowly varying (logarithmic) electrostatic 
potential in two dimensions. 

The solution of the PB equation for the counterion 
density takes different forms depending on whether the 
Manning parameter is larger or smaller than a thresh- 
old value given by A = ln(a out /a)/[l + ln(a ou t/a)] [39i — 
|4~0 | . Here we are interested mainly in the situations with 



sufficiently large Manning parameters £ > A, where the 
normalized counterion density can be written as 



Hp) 



2ir£p 2 



a In 



cot 



(15) 



with a determined from the transcendental equation 



1 



1 — a cot[— a ln(a ou t/o)] 



(16) 



Note that the density profile is renormalized such that 
we have 



2tt 



Hp) p^p = i- 



(17) 



The above density distribution is a monotonically de- 
caying function with very slowly convergent asymptotics 
[50|. It is also independent of the dielectric jump param- 
eter, A, due to the axial symmetry which implies that 
the electrostatic field vanishes inside the cylinder. The 
absence of dielectric discontinuity effects is specific to the 
mean-field limit where no fluctuations are taken into ac- 
count and can be derived also in the case of planar slabs 

S3- 

When the Manning parameter is decreased, the system 
exhibits a continuous counterion condensation transition 
[69| at a critical Manning parameter £ = 1. The nature of 
this transition has been analyzed throughly by means of 
analytical and numerical methods 31( . It was in particu- 
lar shown that the behavior close to the transition point 
is described by the mean-field theory and fluctuation and 
correlation effects (that will be important eventually at 
sufficiently large Manning parameters) play no role at 
the transition point. Here we shall not consider the be- 
havior of counterions close to the transition point but 
do note that since the image-charge effects are absent in 
the mean-field limit, such effects are expected to have no 
influence on the counterion condensation transition itself. 

In what follows, we shall also take into account the 
hard-core repulsion between counterions and the cylinder 
by assuming that the counterions have a finite radius R c . 
Within the PB theory, this can not be done explicitly 
unless full corrections due to excluded- volume repulsions 
between counterions are taken into account ff^L Wty . Here 
we shall consider a simplified version of the ion size effects 
by taking an effective hard-core radius for the cylinder, 
i.e., by setting a —> a + R c . Though approximate on the 
WC level, this procedure turns out to be exact in the SC 
limit (H. 



B. Strong coupling limit 

In the regime when 2 3> 1, the mean-field approxi- 
mation breaks down, and a different kind of approach is 
needed. We shall focus on the limit of S — s- 00, where the 
system can be treated via a systematic strong coupling 
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FIG. 2: a) Radial counterion number density profile around a single charged cylinder as predicted by the SC theory, Eq. (|23l) . 
for various values of the coupling parameter. Vertical dashed line represents the hard-core exclusion volume a + R c , where 
counterions radius in all the figures is taken as R c — 0.2 a. b) Log-linear plot of the counterion density as obtained from the WC 
theory (Eq. (|15[1 . red line), the SC theory (Eq. (|23p . black line) and MC simulations (symbols) in the absence of image-charge 
effects (A = 0). c) Same as (b) but in the presence of a dielectric discontinuity at the cylinder surface with A = 0.95. The 
system is considered in a cylindrical confining cell of radius a out /a = 100. 



virial expansion 0, [I?} leading to an exact analytical 
theory in the leading order in S. By construction the 
SC theory is a single particle theory and contains con- 
tributions only from single-particle interactions. We do 
not elaborate on the derivation of the SC theory here 
as further details can be found in previous publications 

In this Section, we shall be interested in the counterion 
density profile around the cylinder in the SC limit. The 
general form of the counterion distribution within the 
SC theory turns out to be given in terms of a Boltzmann 
factor containing the interaction energy of an individ- 
ual counterion with a fixed charged macroion [14|. In a 



dielectrically inhomogeneous system, such as a charged 
cylindrical core considered here, the SC counterion den- 
sity may be written as [13] 



n(r) cx exp [-/3W se u (r) - /3W 0c (r)} 



where 



-/?(e (?) 2 u im (r,r), 



(18) 



(19) 



is the dielectric image self-energy of a single counterion, 
i.e., the electrostatic energy of a charged point particle in 
the vicinity of a neutral dielectric cylinder. This energy 
of course represents the contribution from the interaction 
of a counterion with its image charge and thus depends 
crucially on the value of the dielectric discontinuity at 
the macroion surface. Furthermore, 



/3W 0c (r) = f3e q / u(r, r')a(r')dr' = u (r) + u im (r), 



(20) 



where 



v Q (r) = (3e q J uo(r > r>(r / )dr / = 2£]np, (21) 

is the bare electrostatic interaction energy of a single 
counterion located at a given position r = (p, tp, z) with 



the cylinder charge (up to an irrelevant additive con- 
stant). It goes to zero when the charge density on the 
macroion surface goes to zero. The image-dependent part 
of the interaction energy between a point charge and 
the surface charge distribution on the dielectric cylinder, 
i>i m (r), can be shown to be nil due to the axial symmetry, 
i.e. 



i'ii 



/3e q 



.(r.r'^rOdr' = 0. 



(22) 



If the surface charge distribution on the macroion is non- 
uniform, such as in the case of charged helical stripes 
[37| . the image dependent part of the interaction energy 
is non-zero and should be considered explicitly. 

Note that physically the dielectric image-dependent 
part of the interaction energy corresponds to the interac- 
tion between the dielectrically induced image charge of 
the surface charge on the cylinder with the counterion as 
well as the image charge of the counterion with the sur- 
face charge on the cylinder, as follows from the general 
definition of the Green's function ([7]). 

Using the above equations, we obtain the SC counte- 
rion number density profile (p > a) as 



n(p) = Cp «exp -rl(p/a 



(23) 



where we have introduced the dimensionlcss integral 

J (x) = E / U{t)K 2 m {tx)dt, (24) 

m=0 J ° 

and the numerical prefactor C is determined from the 
normalization condition, Eq. (|17p . We can extract two 
limiting behaviors for the function I{x), viz. 

7rA 

x -> 1+, 



I(x) 



4(x-l) 

tt 2 A(4 + 3A) 
, 32(1 + A)a; 3 



(25) 



x — > 00. 
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The limiting form of I(x) for x — > 1 + implies that at 
very small separations between a single counterion and 
the cylindrical surface, the dielectric self-energy has the 
same form as in the case of a counterion next to a planar 
wall, where the polarization effects can be described by 
an image charge inside the wall [37| ■ 

The SC density for A = reduces to the well-known 
form for a homogeneous system n(p) oc p~ 2 ^ [3l|. If the 
dielectric core has a smaller dielectric constant than the 
medium (A > 0), as one typically encounters in the case 
of biological macromolecules, the image charges have the 
same sign as counterions and thus lead to a depletion of 
counterions in the vicinity of the cylinder. This behav- 
ior can be seen from the SC density profiles plotted in 
Fig. [2^, as the dielectric mismatch parameter A or the 
coupling parameter S is increased. They show a reduced 
density close to the cylinder (depletion zone) and a peak 
some distance away from the cylinder surface. One can 
estimate the location of the density peak by setting the 
density (|2"5)l derivative to zero, dh(p)/dp = 0. Using the 
approximation (|25p for p = a + h, where h is the peak dis- 
tance from the cylinder surface, we obtain to the leading 
order for small h 



h\ 2 



AS 



(26) 



which shows good agreement with the peak location val- 
ues of the SC density profile as seen in Fig. [2ja,. 

As noted before, effects of a finite counterion radius 
R c within the SC approximation can be taken into ac- 
count exactly (due to the single particle nature of the the- 
ory) by increasing the hard-core radius of the cylinder to 
a — > a + R c . The counterion-counterion excluded- volume 
repulsion effects are absent within the leading-order SC 
theory and enter only in the subleading terms. 



C. Comparison with MC simulations 

Wc have performed extensive MC simulations in order 
to verify the validities of the weak- and strong-coupling 
approaches. The detailed description of MC simulation 
is given in Appendix |Bj 

In Fig. [2p we show the radial density profile of counte- 
rions around the cylinder without dielectric image effect, 
A = 0. Apparently in this case the simulation data are 
always bracketed by the two analytical results obtained 
from the WC and SC theories, Eqs. (flip)) and (j2"5j) . respec- 
tively. The WC theory is found to be valid at sufficiently 
small coupling parameter or sufficiently large radial dis- 
tances from the cylinder. The SC theory is valid in the 
opposite regime, i.e., for sufficiently large coupling pa- 
rameters or sufficiently small radial distances. In fact, 
in agreement with the general trend obtained for planar 
surfaces [II M, 03, the validity regime of the SC (WC) 
theory expands to large (smaller) separations as the cou- 
pling becomes larger (smaller). In the case of a planar 
charged wall, the validity regime of the SC theory can 



be estimated systematically as z/p < VE, where z is 
the distance from the charged plane [lH . In the case of a 
charged cylinder, a similar criterion is yet to be obtained. 

In Fig. [5J:, we show the results for the case where the 
dielectric constant of the cylinder is different from that of 
the medium such that A = 0.95. As already discussed the 
dielectric discontinuities have no effects in the WC limit 
as confirmed also by the agreement obtained between the 
simulation data and the WC (mean-field) result for small 
5. By contrast, the image effects become quite signifi- 
cant at large coupling parameters as the density profile 
deviates qualitatively from those obtained with A = 
(Fig. [2Jd). The simulation data again show good agree- 
ment with the SC approximation at high enough cou- 
plings and small enough radial distances and in particu- 
lar exhibit a depletion zone and a peak at small distances 
from the cylinder surface as predicted within the SC the- 
ory. Note that in a diclcctrically inhomogeneous system 
(A > 0) the SC theory has an explicit dependence on the 
coupling parameter S because of the self-interaction of 
counterions with their image charges. 



III. TWO LIKE-CHARGED CYLINDERS 

We now turn our attention to the problem of the in- 
teraction between two like-charged cylinders, which is 
commonly used as a model for the interaction between 
rigid polyelectrolyte chains [38[. We consider two identi- 
cal parallel cylinders of radius a oriented along the z axis 
at an interaxial separation of R (see Fig. |3|). The elec- 
tric charge is assumed to be uniformly distributed on the 
surface for both cylinders with equal linear charge den- 
sity A. The cylinders are assumed to be infinitely long 
of length L — ¥ oo and are standardly confined within a 
confining square box of lateral size L± (we shall assume 
L±/a = 60 in Section fill Al and L±/a = 100 in Sections 
IIII BllrnD|) . but we should emphasize that in the regime 
of Manning parameters considered here, the lateral box 
size plays no significant role [3l|, HfjJ . 

We follow the same approach as in the case of a single 
cylinder and thus first analyze the problem in the weak- 
coupling limit and then in the strong-coupling limit. We 
then compare the results from these two analytical theo- 
ries with the MC simulations. 



A. Weak coupling limit: Mean-field theory 

The weak coupling regime (3 <C 1) is again charac- 
terized by the mean-field PB equation. In this case how- 
ever, no closed- form solution is available for the nonlinear 
PB theory. We thus take full recourse to the numerical 
methods appropriate to the problem. In the two-cylinder 
geometry, the axial symmetry is broken but the transla- 
tional symmetry along the z axis remains intact; hence 
one has to deal with a two-dimensional PB equation in a 
finite bounding square Lj_ x L±, which is then solved nu- 
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FIG. 3: (a) Schematic top view of two identical charged di- 
electric cylinders oriented parallel in the z direction. The 
system is confined in a square box of lateral size L± contain- 
ing the two cylinders and their neutralizing counterions (not 
shown), (b) The potential generated by two charged dielec- 
tric cylinders can be reproduced equivalently by four linear 
(image) charges as explained in the text. 



merically [74J . Wc shall focus on the interaction force be- 
tween the cylinders, which can be evaluated via the elec- 
trostatic stress tensor in a standard manner [6(|. Note 
that the WC scheme used here based on the numerical 
solution of the PB equation allows for a full analysis of 
the dielectric discontinuity effects. (Later on in our SC 
analysis, which will be based on analytical methods, we 
shall restrict our discussion to only the first-order image 
approximation, as a full analysis of dielectric discontinu- 
ity effects in the SC limit is not yet available.) Here again 
the counterion size effect will be taken into account ap- 
proximately by setting the distance of closest approach 
to each cylinder as a — > a + R c . Fluctuation contribution 
around the mean-field solution in cylindrical geometry is 
non-trivial to calculate [49] but is by construction small 
and we thus skip its detailed analysis. It does however 
depend crucially on the dielectric discontinuity at the 
cylinder surface. 

In general, the WC force obtained between two like- 
charged cylinders is repulsive and decays monotonically 
with their interaxial separation R. The mean-field pre- 
diction for the rescaled force per unit length defined via 



F = /2££n\ F 



(27) 



is shown in Fig. [4] as a function of the interaxial separa- 
tion R and for different values of the dielectric disconti- 
nuity parameter A (see Eq. @). Note that in contrast 
to the case of a single cylinder studied in the previous 
Section, the dielectric discontinuity is expected to mat- 
ter in the two-cylinder geometry even within the mean- 
field approximation. This is due to the axial symmetry 
breaking which leads to electric field penetrating into the 
cylindrical cores. As seen in the figure, the WC mean- 
field results compare very well with the MC data (to 
be discussed later) at sufficiently large interaxial separa- 
tions. The dielectric discontinuity has a relatively small 




FIG. 4: Rescaled force between two like-charged dielectric 
cylinders as a function of their interaxial separation. The WC 
prediction obtained by solving the corresponding PB equation 
(solid lines for A = and 0.95) are compared with the simu- 
lation data (symbols). 



effect on the interaction force. It leads to an increased 
repulsion (and also larger deviations from the mean-field 
prediction) at small interaxial separations but its effects 
diminish at larger separations. 



B. Strong coupling limit 

As noted before, the SC theory follows systematically 
from the leading order contribution to the partition func- 
tion in the limit of infinite coupling parameter S — > oo 
fl4| . The corresponding SC free energy involves only 
interactions between single counterions and the charged 
cylinders as well as the direct electrostatic interaction 
between the cylinders themselves (see Refs. [l4|, EH). In 
the presence of the dielectric discontinuity effects, the SC 
free energy per counterion can be written as 



N 



N 



In / exp [-/3W se if (r) - /3W 0c (r)] dr. 



(28) 

Here N is the total number of counterions and Woo is the 
electrostatic energy due to the interaction between the 
two cylinders in the absence of any counterions, which 
contains both the direct Coulomb interaction between 
their surface charges as well as the contribution due to 
their images that results from the polarization of their 
dielectric cores. Furthermore, Wq c is the energy due to 
the interaction between a single counterion and the sur- 
face charges on both cylinders and includes contributions 
from the corresponding image charges as well. Finally, 
Wgcif is the image self-energy of counterions, i.e., the 
contribution from the interaction of an individual counte- 
rion with its own image charges in both cylindrical cores. 
Note that the volume integration should be again taken 
over the space available to counterions, i.e., inside the 
confining square box of lateral size L± excluding the two 
cylindrical cores. Note also that on the SC level, the 
counterion excluded- volume repulsions are absent (as im- 
plied by the fact that strongly coupled counterions are 
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highly isolated within large correlation holes in the SC 
limit) and only the excluded-volume interaction between 
the counterions and the cylinders will be present. In fact, 
the counterion size effects can be accounted for exactly 
within the SC theory via a hard-core repulsion which sim- 
ply amounts to setting the effective (hard core) cylinder 
radius as a + R c [6(| ■ 

As in the single cylinder case all charges are coupled 
with the interaction kernel, composed of direct and im- 
age part, but with the difference that now the image 

(2) 

part of the kernel W; m ; (r,r') corresponds to two dielec- 
tric cylinders separated by the distance R. This kernel 
should satisfy the electrostatic boundary conditions on 
both cylindrical surfaces which leads to very com- 
plicated numerical expressions and is unfortunately not 
available in a closed analytical form. Our aim is never- 
theless to give a simple approximate analytical expression 
for the final results, so we need to make an approximation 
at this step. In what follows we focus on the first-order- 
image approximation, where we neglect higher orders of 
inter-cylindrical image interaction. This approximation 
implies that the image Green's functions of the two cylin- 
ders can be written as the sum of the Green's functions 
of the isolated cylinders, Eq. 



u im( r > r ') ~ «im(ri,ri) + u im (r 2 ,r' 2 ), 



(29) 



where Ui m (r,r') is the one-cylinder image kernel, 
Eq. (fTTj) . Here, ri, r[ are distances centered at the first 
and r 2 , r 2 at the second cylinder, for instance, 



n, 2 = r T (R/2)e x , 



(30) 



with e x a unit vector pointing in the x direction, i.e., 
horizontally on Fig. [3] Recall that the surface charge is 
uniformly distributed on both cylinders so that the fixed 
charge distribution is composed of two single cylinder 
distributions 



(31) 



where cr(r) is the charge distribution of a single cylinder, 
Eq. ([5]). With these assumptions one can then evaluate 
all interaction terms in Eq. (f2"8"| analytically as follows. 

The cylinder- cylinder interaction energy per counte- 
rion can be obtained as 



^T = 5^ //^)(ry 2 )(r,r')^)(r')drdr', 



(32) 



Furthermore the image self-energy of a single counterion 
interacting with the dielectric cores of both cylinders can 
be derived in the form 

/3W se if (r) = i/3(e (7) 2 ii im (ri, ri) + i/3(e <7) 2 Wi m (r 2 , r 2 ) 



[I( Pl /a)+I( P2 /a)} 



(33) 



where I(x) is defined in Eq. (|24|) . Finally, the single par- 
ticle cylinder-counterion contribution can be written as 
the sum of the interaction with each cylinder separately 
and assumes the form 

/3W c(r) = Pe q J u^(r, r') a< 2 >(r')dr', 

= w(ri) + w(r 2 ). (34) 

The cylinder-counterion energy contribution from the fx- 
th cylinder (with /i = 1,2) is itself composed of three 
parts, 

u(r M ) = Uo(l>) + Wsamc(r Al ) + Wcross(r A1 ), (35) 

which correspond respectively to three parts of the 
Green's function, i.e., u^ 2 ^(r,r') « Uo(r, r')+u; m (ri, r' 1 ) + 
Wim(r 2 , r 2 ), within the approximation implied by Eq. (|2"9")) 
and will be evaluated as follows. 

The first contribution in Eq. (|3"5j) , wo(r M ), is the direct 
interaction of the bare surface charge of each cylinder 
with a single counterion, i.e., for the the p-th cylinder, 

«o(i>) = (3e q /"^(r^r^o-fr^dr', 

= 2£]np„, (36) 

where p^ is the radial distance from the center of the 
/i-th cylinder to the counterion, and £ is the Manning 
parameter as defined in Eq. (U) . The second contribution, 

Usame(r M ) = /3e a q J u im (r M , rj i )<7(rj 4 )dr' = 0, (37) 

corresponds to the interaction of a counterion image 
charge with the surface charge on the same cylinder, 
which is thus zero by symmetry reasons as already dis- 
cussed in the one-cylinder case, Eq. (f2"2"j). This contri- 
bution can be also considered as the interaction between 
the image of the surface charge on the same cylinder with 
the counterion. 

Finally, we have the cross contributions due to the in- 
teraction of a counterion image in one cylinder with the 
surface charge on the other cylinder. Formally, it follows 
for each cylinder as (see Appendix fA"f 



Ucross(i>) =0e o q J u im (r AI ,r^)cr(r' l/ )dr', 

= 2A£ln(^Y (38) 



for /j/v, where is the radial distance from an axis 
shifted by a 2 /R from the center of the /x-th cylinder to- 
ward the axis of the other cylinder, sec Fig. [3b, i.e., 



, 2 x2 



Pi =P,-2(-jp„cos^ , 

or in cartesian coordinates, 

-i? _ a 2 
2 ~ ~R 



±x) +y. 



(39) 



(40) 
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FIG. 5: Rescaled SC force between two like-charged cylinders 
in the absence of a dielectric discontinuity (A = 0) as a func- 
tion of the rescaled interaxial separation for various Manning 
parameters as indicated on the graph. 



The SC free energy thus follows by inserting the above 
expressions into Eq. (|28[) in the form 



N ^ ^ \ R 2 J 



-In 



■£(1 - A)lnp?p2 _ £AlnpJ 2 p| 2 dr. (41) 



The corresponding SC force between the two cylinders is 
obtained simply by differentiating the free energy with 
respect to R 



F = - 



dr_ 

dR' 



(42) 



which may be expressed in dimcnsionless units according 
to Eq. (pf|) . 

The SC density profile follows from the general expres- 
sion (|18p . which is the same as the integrand in (|4Tj) . viz. 



n(r) = Ccxp<^ 



(43) 



£(l-A)lnp 2 p 2 -£Aln Pl 2 ^ 



where the normalization prefactor is determined from the 
electroneutrality condition. 

Note that the forms of the SC free energy, Eq. (|4*T|). 
and the SC counterion density profile, Eq. (|43j) . reflect 
the fact that, within the first-order image approximation 
(and using the appropriate volume constraints for coun- 
tcrions), the two charged dielectric cylindrical cores can 
be equivalcntly replaced by four parallel lines of charges, 
two of which are located along the central axis for each 
of the cylinders with (renormalized) linear charge density 
(1 — A) A and the two other lines of charges with linear 
charge density +AA are placed along two axes shifted by 
a 2 /R from the central axis for each of the cylinders, see 
Fig. |31 Note that the image charge language is applica- 
ble for two charged dielectric cylinders in the absence of 
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FIG. 6: Counterion density around two parallel like-charged 
cylinders with Manning parameter £ = 10 at interaxial sepa- 
rtaion R — 3a and in the absence of image charges (A = 0). 
The density is shown across an arbitrary plane perpendicular 
to the two cylinders in a color-coded fashion as indicated on 
the graph. The small depletion zone around each cylinder 
is due to the counterion volume exclusion as counterions are 
assumed to have a finite radius of R c — 0.2a. 



pointlike charges (counterions) [701, an d can be exactly 
solved by applying higher orders of images. But in order 
to be consistent with the first-order image counterion- 
cylindcr approximation, we have used the same image 
kernel (|29[) . which accounts for the first-order image ap- 
proximation between the two cylinders as well 



C. 



SC interaction in a dielectrically homogeneous 
system (A = 0) 



In order to proceed, let us first consider the case where 
there is no dielectric mismatch between the cylinders and 
the surrounding medium, i.e., A = 0; the image charges 
will thus be absent. The SC interaction in this case has 
been considered also in previous works [3(| IH, H3] . We 
shall reproduce some of these previous results for the 
sake of completeness, but will also provide several new 
results, including a direct comparison with MC simula- 
tions for the force dependence as a function of separation 
at different values of the coupling parameter as well as 
a global interaction phase diagram, which have not been 
considered previously. 

In the absence of image charges, the SC free energy 
(|4"Tj) reduces to a simple form as 



— = -£mi?-ln 
N ^ 



J exp (-2£_\np 1 p 2 ) dr, 



(44) 



where the first term is again the bare repulsion between 
the two cylinders and the second term contains energetic 
and entropic contributions from counterions. The results 
for the SC force that follow from Eq. (|44p are shown in 
Fig. [5] for a few different Manning parameters, where we 
have accounted for the finite counterion radius by choos- 
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FIG. 7: Rescaled interaction force between two like-charged cylinders in the absence of image charges (A = 0) as a function 
of the rescaled interaxial separation for Manning parameters £ = 2 (a), 10 (b), 100 (c) and different values of the coupling 
parameters as indicated on the graphs. The symbols represent the simulation data, the red lines the prediction of the SC theory 
(Eq. (|44)0 and the blue lines are that of the WC theory ( Section UlI A[) . 



ing R c = 0.2a. The SC force is repulsive at small sep- 
arations and becomes attractive beyond the equilibrium 
interaxial distance R* where the force vanishes. The SC 
attraction is mediated by countcrions that are strongly 
coupled to both cylinders in this limit and are thus accu- 
mulated mainly in the region between the two cylinders 
as can be seen from the counterion distribution (|43|) in 
Fig. |51 This is a direct consequence of the energy contri- 
butions included in the second term in Eq. (|44p j30L H|| . 

By decreasing the separation between the cylinders, 
the entropic osmotic pressure from counterions sand- 
wiched between the two cylinders becomes increasingly 
important and the effective interaction becomes repul- 
sive. This repulsion is reduced as the surface-surface 
separation, R — 2a, becomes smaller than the counterion 
diameter, 2R C as the counterions are depleted from the 
intervening region due to excluded-volume effects. This 
behavior due to the presence of the counterion depletion 
interaction has been investigated thoroughly in the case 
of cylinders with helical charge pattern and we shall not 
delve on it any further here [60| . 

As seen in Fig.[5]thc entropic repulsion effects at small 
separations decrease when the Manning parameter £ is 
increased. For very large Manning parameters, the free 
energy is dominated by purely energetic contributions 
and does not contain any temperature effects. In fact 
the limit £ — > oo formally corresponds to the zero tem- 
perature limit as we have already taken the limit 3 — > oo 
appropriate within the SC approximation. The connec- 
tion between the SC limit and the zero temperature limit 
has been analyzed in detail in Refs. [HI, 53- The asymp- 
totic form of the SC force in the limit of large £ can be 
obtained as 



F/L 



4tt£ 
~Rj~a~ 



47r£/_ 
R/a\R 



2R 



Rc 



R<2(a + R c 



1 R>2(a + R c ), 



librium hard-core surface-surface separation 
SR* = R* -2(a + R c ), 



(46) 



tends to zero as SR*/a ~ 2/(30 as £ -> oo gfj. Note 
that even for moderate values of the Manning param- 
eters, £ ~ 5, the SC theory predicts a closely packed 
bound state with a relatively small surface-surface sepa- 
ration SR* I a ~ 0.1. 

Let us now consider the situation where the coupling 
parameter has a finite value. In this case we shall study 
the system using MC simulations as detailed in Appen- 
dices [B] and The MC results for the interaction force 
between like-charged cylinders (with no dielectric mis- 
match) are shown in Fig. [7] along with the WC and SC 
predictions. As seen in Fig. [7^,, the WC results agree 
with the simulation data for a sufficiently small coupling 
parameters 3 < 1 as expected. Furthermore, by increas- 
ing the coupling parameter, the MC results start deviat- 
ing strongly from the WC prediction and tend to the SC 
prediction, where a reasonable agreement is obtained al- 
ready at coupling parameters of the order 5 ~ 100. Note 
that the WC and SC predictions bracket the simulation 
data and thus establish upper and lower bounds for the 
interaction force at any realistic value of the coupling 
parameter. 

In Figs. [TJd and c we show the results for larger Man- 
ning parameters of £ = 10 and 100, respectively. As Man- 
ning parameter becomes larger, a larger 3 is required in 
order to achieve the same level of agreement with the 
SC prediction at a given separation and thus the con- 
vergence to the SC limit becomes weaker. This effect 
is more pronounced at larger interaxial separations. In 
fact, as is generally known [14|, the SC effects become 



(45) 

which is shown as a black solid line in Fig. [SJ The equi- 



more dominant and the SC theory becomes increasingly 
more accurate at smaller surface-surface separations. In 
general, the applicability regime of the SC theory (ob- 
tained in the limit of 3 — > oo) to the situations where 
the coupling parameter 5 is finite can be specified via 
simple validity criteria as applied successfully in several 
previous studies EE H IS IM El E3 • In what fol - 
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(a) (b) (c) 

FIG. 8: (a) Schematic representation of the validity criterion for the SC theory at large £. The SC predicts the configuration 
where all counterions are located between the cylinders. This argument fails if counterion-counterion repulsion is too strong, 
(b) Bound-state interaxial separation as a function of the Manning parameter. The solid line shows the SC prediction and the 
symbols represent the simulation data. The horizontal dashed line represents the closest approach distance between the two 
cylinders, 2a, and the vertical dashed lines represents the minimum value of the Manning parameter = 2/3 where SC force 
may be attractive, (c) Phase diagram indicating the repulsion and attraction between two like- cylinders with no image charges 
in terms of H and £ in the absence (blue squares) and in the presence of image charges (red triangles with A = 0.95). Symbols 
show the boundary line obtained from MC simulations. Filled circles show the points appropriate for DNA with counterions 
of valencies in the range q = 1, . . . , 6. In (b) and (c) dashed lines are guides to the eye. 



lows wc shall derive the particular form of the SC validity 
criterion for the case of two like-charged cylinders. 

It should be noted first that the SC theory (obtained as 
an exact asymptotic theory in the limit of S — ¥ oo) con- 
tains only contributions stemming from the interaction 
of individual counterions with charged objects and thus 
the counterion-counterion repulsions and other higher or- 
der many-body effects are absent. These effects however 
become increasingly more important as the coupling pa- 
rameter is decreased, depending crucially on the precise 
value of all the system parameters. In the case of two 
cylinders, as noted before, counterions are mostly accu- 
mulated in the intervening region between the two cylin- 
ders (Fig. |6|) where the typical counterion spacing can be 
estimated by stipulating the local electroneutrality con- 
dition as d = qeo/(2X) or in rescaled units as d/(i = S/2£. 
It is thus evident that as the Manning parameter becomes 
larger at a fixed coupling parameter the spacing between 
counterions tends to become small. Hence, a larger effect 
due to counterion-counterion repulsions and thus larger 
deviations from the SC theory would be expected to be 
observed at larger £. In other words, the SC theory is 
expected to overestimate the counterion density in the 
intervening region at any finite value of the coupling pa- 
rameter, which is consistent with the general result that 
the SC theory gives the upper bound for the density pro- 
file of the counterions [HI, [ljl [l?], HI] . One can thus iden- 
tify the validity regime of the SC theory by estimating 
the counterion-counterion contributions and comparing 
them with the sinlge-particle counterion-cylinder contri- 
bution at a given interaxial separation. This can be done 
by noting that the ground state configuration (obtained 
for £ — > oo and S — > oo) is predicted within the SC theory 
to be the configuration where all counterions are localized 



between the two cylinders and are lined-up along the z 
axis as obviously favored by the electrostatic interaction 
energy. 

When the coupling parameter is reduced to a finite 
value, repulsions between counterions become more im- 
portant and some counterions tend to escape from the 
intervening region to the external side of the cylinders 
(shown by an arrow in Fig. |8ji). In order to prevent 
this from happening, the coupling parameter 3 has to be 
large enough. By comparing the energy of the ground 
state with the energy of a single defect obtained by al- 
lowing one counterion to move to the external side of one 
of the cylinders, it follows analytically that for a large £ 
and in the absence of image charges, one should have 



E > 4.2 £ 2 



R ~ 2a, 



S>2.9£ 2 (^) i?>2a. 



(47) 



These criteria indeed capture the trends observed in 
Fig. [7] when viewed in terms of the three main param- 
eters describing the system, i.e., 2, £ and R. They au- 
tomatically cover the one proposed before based on the 
surface-surface separation between the cylinders, where 
one requires that the surface-surface separation between 
the cylinders must be smaller than the s pac ing between 
counterions, SR < d, or in rescaled units |30j, |46[ 



5R/[i < S/(2£). 



(48) 



This latter criterion is inspired by the observations in 
the case of charged planar walls where it can indeed be 
derived from systematic analysis of higher-order correc- 
tions to the SC theory 14|. It can describe the validity 
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regime of the SC predictions for the equilibrium sepa- 
ration between strongly coupled cylinders that form a 
closely packed bound state [30, 56j, but when the ap- 
plicability of the SC theory for the interaction force is 
considered, the criteria (|47|) should be used instead. The 
difference between the new criteria (|47[) and the one in 
Eq. (UHJ) is in fact related to the fundamental difference 
between the topology of the space available to counte- 
rions in the case of two cylinders as compared to two 
planar walls: in the latter case counterions can not "es- 
cape" from the intersurface gap and the ground state of 
the system has a universal two-dimensional configuration 
HH EE G3 , while in the former the ground state corre- 
sponds to a one dimensional arrangement of counterions 
and the whole space is available for thermal excitations 
from the ground state. 

The criteria (|47[) are more stringent and cover also the 
situation where the cylinders are placed at large interax- 
ial separations. The ground state of neutralizing coun- 
terions between two charged cylinders (without the ef- 
fects of the dielectric mismatch) has been investigated 
extensively by Arnold and Holm in Ref. [47j ], where a 
validity criterion for the SC theory has been proposed 
as 3 > 3.45 £ 2 (i?/a) (when expressed in rescaled units) 
based on computer simulations. This agrees with our es- 
timate in Eq. (|4T|) within 16% of the numerical prefactor. 

Let us now turn our attention to the behavior of the 
equilibrium interaxial separation or the so-called bound- 
state separation, R* , where the interaction force between 
the two cylinders vanishes. As noted before, the SC the- 
ory is generally more accurate at smaller separations. A 
closer inspection of Fig. [7] shows that the SC prediction 
for the bound-state separation agrees with the MC data 
even outside the regime set by the SC criteria (|47)) , which 
reiterates the point discussed above that a less stringent 
criterion such as Eq. (|4"5)l would be sufficient to describe 
the validity of the SC prediction for the bound-state sep- 
aration. In Fig. [8]d the bound-state interaxial separation 
is plotted as a function of the Manning parameter, which 
shows that a reasonable agreement between MC data and 
the SC theory for this quantity can be achieved for a cou- 
pling parameter as small as 2 ~ 50, in agreement with 
previous results in Ref. (30|. As seen in the figure, the 
bound-state separation diverges (i.e., the cylinders un- 
bind) for Manning parameter approaching a minimum 
value of £* = 2/3 [46[ below which the two cylinders 
merely repel each other. 

Figure [8)3 can be also thought of as a phase diagram 
identifying the attraction (the region above the lines) and 
repulsion (the region below the lines) regimes between 
two like-charged cylinders in terms of the parameters R* 
and £ for a given coupling parameter 2. 

In Fig. [5J: we show a global phase diagram in terms 
of the parameters 2 and £. The boundary lines are de- 
termined from MC simulations (open symbols) and sep- 
arate the region (above the lines) where the interaction 
force becomes attractive at some finite interaxial sepa- 
ration between the cylinders and the region (under the 



(c) 5 = 1000 (d) 5 = 5000 

FIG. 9: Counterion density around two parallel like-charged 
cylinders with Manning parameter £ = 10 and the dielec- 
tric discontinuity parameter A = 0.95 at interaxial separ- 
taion R = 3a. The density is shown across an arbitrary plane 
perpendicular to the two cylinders in a color-coded fashion 
as indicated on the graph. The small depletion zone around 
each cylinder is due to the counterion volume exclusion as 
counterions are assumed to have a finite radius of R c = 0.2a. 
A large depletion zone develops due to the influence of image 
charges as the coupling parameter increases from H = 100 to 
5000 (a-d). 



lines) where the force never becomes attractive at any 
interaxial separation. In other words, the boundary lines 
themselves correspond to the parameter values where 
the force-distance diagrams exhibit a global minimum 
at F = 0, i.e. the force curve only touches the abscissa. 
They were determined by bisection procedure for a series 
of MC runs. If the cylinder parameter values are chosen 
to describe DNA, we have £ = 4.1 q and 2 = 2.8 q 3 for 
q valency counterions, which are shown as solid circles 
in Fig. \8jp. This suggests that DNA-DNA interaction in 
the case of monovalent counterions falls well within the 
repulsive region, while with divalent counterions and be- 
yond one can expect to observe an attractive force, in 
qualitative agreement with recent experiments @~E1 ■ 



D. Image-charge effects (A > 0) in the SC limit 

We now consider the case of two charged cylinders with 
a dielectric core whose dielectric constant may be in gen- 
eral different from that of the surrounding medium. In 
this case, the induced image charges (due to the polariza- 
tion of the dielectric cores) can play an important role as 
well. However, as we discussed before (Fig.@|, the image- 
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(a) (b) (c) 

FIG. 10: (a) and (b) Force results between two dielectric charged cylinders (A = 0.95) for various coupling parameters S. The 
MC results (symbols) are compared with the corresponding (same color) SC results obtained from Eq. (|41|l and shown here as 



solid lines. In (a) the MC results for S = 1 are compared with WC prediction (solid green line) as discussed in Section [III Al 
The dashed lines shows the bare electrostatic interaction between the two cylinders in the absence of counterions, Eq. (|C3|I in 
the appendix, (c) Bound-state distance for two cylinders with A = 0.95. Here MC results (symbols) are compared with the 
SC predictions (solid lines of the same color) . 



charge effects turn out to be small in the WC regime 
where the coupling parameter is small (e.g., when coun- 
terions arc monovalent). This may be the reason why 
such dielectric effects have not yet been throughly inves- 
tigated in the particular case of two charged cylinders. 
The situation turns out to be quite different in the oppo- 
site limit of strong coupling as we will see in this Section. 
This is tied in with a fundamental difference between 
countcrion distribution in the SC limit as opposed to the 
WC limit: while in the former limit one deals with highly 
isolated single counterion close to a charged surface, in 
the latter limit individual nature of counterions fades in 
the wake of dominant collective many-body effects (as 
counterions form a diffuse and uncorrelated ionic cloud 
around the cylinders) and thus the polarization effects 
are highly reduced. 

If the dielectric cores have a smaller dielectric constant 
than the medium (A > 0), as relevant to most macro- 
molecules in water, the induced image charges will have 
the same sign as the counterions and thus tend to cause 
depletion of counterions from the vicinity of the dielec- 
tric cores. This behavior is shown in Fig. El where the SC 
density is plotted across an arbitrary plane perpendicu- 
lar to the two cylinders in a color-coded fashion. As can 
be seen, the depletion effect is enhanced as the coupling 
parameter, S, is increased such that a large "depletion 
zone" develops that eventually encloses both cylinders. 
(Note that the SC free energy in the presence of a dielec- 
tric discontinuity, Eq. (|41[) . has an explicit dependence 
on the coupling parameter E as alread y k nown from the 
case of planar charged dielectric slabs [37[). 

In fact, the image depletion at larger couplings turns 
out to have a similar effect qualitatively as the counte- 
rion volume exclusion, i.e., as if the counterion radius R c 
is renormalized to a large effective value. Note that the 
region with highest concentration (shown in red color) is 
"squeezed" and eventually splits in two regions shifted 
away from the common x — z plane that passes through 



the axes of the two cylinders. The splitting appears when 
the second derivative of the density (f4"3"]) with respect 
to the y coordinate becomes positive at the midpoint, 
d 2 n(r)/dy 2 \ y= o > 0. Using the small-distance approxi- 
mation p — > a + from Eq. (|25[) . we can estimate the dis- 
tance R at which the density splitting happens as 



R-2a\ 2 



AS 



2(1 + SA)? 



(49) 



Thus for the Manning parameter £ = 10 and at sep- 
aration R = 3a, the splitting is estimated to occur at 
S ~ 800 as is the case in Fig. [9l 

In Figs. ITOk and b, we show the SC force per unit 
length (solid lines) as obtained from Eqs. PX| and P2")l 
along with the MC simulation results (symbols) for two 
different Manning parameters and several different values 
of the coupling parameter. As seen upon increasing the 
coupling parameter, the repulsive peak at small separa- 
tions shifts towards larger values of the interaxial separa- 
tion R, which reflects the emergence of a depiction zone 
due to repulsions between counterions and their image 
charges as discussed before. As a result, the main con- 
tribution to the inter-cylinder interaction at small sepa- 
rations comes from the bare interaction between the two 
cylinders and the force is approximately given by the bare 
repulsion force i*oo as given by Eq. (|C3[) in the appendix 
(dashed line, Fig. [TOb). At larger separations, we still 
find strong SC-type attraction even in the presence of 
image charges. One should thus note that in general the 
repulsive forces are enhanced in the effective interaction 
in the SC limit especially at small separations as com- 
pared to a homogeneous system with no image charges. 
Therefore, unlike in the WC limit (see Fig. [2}, the intro- 
duction of image charges may lead to qualitative changes 
in the behavior for a strongly coupled system. 

In general, the SC results for the interaction force show 
better agreement with the simulations for larger 3 and 
smaller separation distances R. Especially the results 
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for the bound-state interaxial separation show very good 
agreement between the SC results (solid lines) and the 
MC data (symbols) as shown in Fig. ITOb (compare with 

Fig. Eh). 

In plan-parallel geometry the repulsive self-image in- 
teraction compresses counterions toward the mid-plane 
and hence expands the range of the SC regime [H| . Here 
in the cylindrical geometry such mechanism cannot work, 
since counterions can escape from intcr-cylindrical re- 
gion. Therefore, it seems that dielectric image effects 
do not have a noticable influence on the range of validity 
of the SC approximation. 

Note that just as in the case with A = (blue line 
in Fig. [5b ) , the simulation results show no attraction be- 
tween the cylinders at too low electrostatic coupling pa- 
rameters. The attraction can only appear once the cou- 
pling parameter exceeds a threshold value. This thresh- 
old value is shown as a red line in the global phase di- 
agram in Fig. [8b for the dielectric discontinuity param- 
eter A = 0.95. Note that bare cylinder-cylinder repul- 
sion is enhanced in the case of images (A > 0), eq.([28|). 
Therefore, larger coupling parameter S is needed to get 
attraction at smaller Manning parameters compared to 
no-image case (red line stays above blue for £ < 3 and 
A = 0.95 in Fig. Et) - For larger Manning parameters, the 
attractive interaction between the counterion image and 
the cylindrical charges (represented by the term v cr0 ssj 
Eq. (|38p) becomes more important and, hence, attrac- 
tion can appear even at lower coupling parameters (red 
line is below the blue line for £ > 3 on Fig. [5b). 



IV. CONCLUSIONS 

We have analyzed the electrostatic interaction between 
two like-charged cylindrical macromolecules surrounded 
by counterions. We have in particular examined the role 
of image charges when the cylinders have a dielectric 
constant which is different from that of the surrounding 
medium within the the weak- and strong-couplin g fra me- 
works elaborated first by Netz and coworkers [Til Il6l. [l7j . 
The two limits are defined in terms of a single coupling 
parameter S. 

The weak coupling limit, or equivalently the mean- 
field theory, is relevant when the valency of counterions 
and/or the macromolecule surface charge density is small 
and/or when the dielectric constant of the medium is 
large enough. It is based on the Poisson-Boltzmann equa- 
tion which turns out to give purely repulsive interactions 
between two-like charged cylinders. The image-charge ef- 
fects in this case turn out to be small even at small sep- 
arations and quickly diminish as the inter-cylinder sepa- 
ration is increased. 

On the contrary, the strong coupling limit is rele- 
vant when the valency of counterions is large and/or 
the macromolecule surface charge density is large and / or 
when the dielectric constant of the medium is small 
enough. The SC approach is based on a single-particle 



description which is obtained as an exact limiting result 
for large coupling parameters 3 — > oo 0, [l|| [l?} • While 
some aspects of the SC theory for cylindrical macro- 
molecules were studied in previous works [30|, HH, HH, |47| , 
other aspects have remained unadressed. Our work is 
aimed specifically at addressing the effects due to image 
charges, which indeed turn out to be quite significant in 
the SC limit. Using a generalized SC theory and exten- 
sive MC simulations we have shown that the counterions 
are strongly depleted away from the cylindrical cores due 
to repulsion from image charges leading to stronger ef- 
fective repulsions between the two cylinders at smaller 
separations. The counterion-mediated attraction will be 
present at intermediate to large separation when the cou- 
pling parameter is sufficiently large. These features are 
in marked contrast with those obtained with the WC 
theory. We have mapped a global phase diagram where 
attractive and repulsive interactions emerge between two 
like-charged cylinders with or without the image-charge 
effects. 

Here we have employed a first-order image approxima- 
tion to deal with the image-charge effects within the SC 
theory (note, however, that the numerical scheme used 
within the WC limit allows to account for the full effect 
of the dielectric discontinuity) . Generally, the treatment 
of the image charges in nontrivial geometries becomes 
very complicated and the two-cylinder model is no ex- 
ception. The first-order image approximation ([^1 has 
the obvious advantage that it allows for an analytical 
treatment of the system. The study of the full effect due 
to the dielectric discontinuity in the two-cylinder model 
requires more advanced numerical and analytical devel- 
opments that might become available in the future. An- 
other interesting effect which can be investigated in the 
present context is that of the additional salt that may 
be present besides the neutralizing counterions. Recent 
generalizations |36j should allow for a systematic study 
of salt screening effects within the WC-SC framework. 
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Appendix A: Cross interaction term t> C ross(r) 

In this Section, we explain how the cross energy term 
(|38[1 may be evaluated in the case of two charged cylin- 
ders (Scction lHI Bp . This requires employing some math- 
ematical identities as we shall explain below. Note that 
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here the position vectors centered at two different cylin- 
ders appear simultaneously, p ^ v, 



Across V 1 - /J 



(r M ) = /3e g 



r^rXOdr'. 



(Al) 



Inserting the image kernel Mi m (r,r'), Eq. (jlip . and the 
cylinder charge density <x(r), Eq. ([5]), we first transform 
the coordinates p^ and </? M centered at the p-th cylinder 
to the those centered at the f-th cylinder using the Graf's 
addition theorems (75| . 

oo 

sin TYMp^ — ^ ^ S 7nn {kpy) sin rupv, (A2) 

n=0 

oo 

cos n</3i/, (A3) 

with m and n being integers and 

{kp v ) = [K m+n {kR) - K m - n {kR)]I n (kp„) (A4) 
C mn (kp v ) = [K 

m-\-n 

{kR) + K m - n {kR)]I n {k Pv ) 



1 



x (1 - -<W 



(A5) 



We can then carry out the integration in (|A1[) over p' 

and z' . Integration over z' produces a {2smkz ao )/k 
term, where z x — > oo is the upper limit and thus we 
obtain 



oo „oo 
m=0 - 70 

x - — — — C m o(ka) cosrrupp. (A6) 

In the next step we integrate over the wave- vector k, 
where the integrand contains the rapidly oscillating fac- 
tor sinfczoo, which supresses contributions of the inte- 
grand in whole the range, except at k — > 0. Formally, we 
use the following mathematical identity 

/>oo 

hm / /(;)sinw tdi= jRes(/,0), (A7) 

where Res(/, 0) is the residue of function f(t) at t = 0. 
This integral is finite only if f(t) behaves as 1/i for t 
and goes to for t — > oo. This statement can be proven 
by considering the split of function f(t) into a 1/t term 
and a remainder, i.e., f(t) = t _1 Res(/, 0)+/(i). The first 
part can be integrated straightforwardly, whearas the 
remainder yields according to the Ricmann-Lcbcsgue 
lemma [76j]. This finally leads us to 



Across (j-fl ) 



oo 
* — ' m 



cos mtpp, 



(A8) 



m=l ■ 

which may be simplified further using the identity 

,2 

cos mt^ = (A9) 



OO -. 



-In 



a 



R P , 



(- 

\Rp, 



and yields the result given in Eq. (|38 



Appendix B: MC simulations 

We performed Monte-Carlo simulations in order to 
study the system of one and two charged dielectric cylin- 
ders at arbitrary coupling parameters, which thus go be- 
yond the weak and strong coupling theories. All simulta- 
tions were performed in the canonical ensemble (NVT) 
using standard Metropolis algorithm [jj [zll • in the case 
of one charged cylinder the system (including the neutral- 
izing counterions) is enclosed in a cylindrical simulation 
box of outer radius a out = 50a whereas in the case of 
two cylinders the system is enclosed in a square box of 
lateral size L^/a = 60 (WC) and 100 (SC) (note that in 
this work we focus on the regime of large Manning pa- 
rameters where the lateral size and shape of the confining 
box becomes unimportant [3ll |4(| ) . The simulation box 
height is assumed to have a finite value L z , which in terms 
of other physical parameters by invoking the global elec- 
troneutrality condition: in the case of a single cylinder it 
is given by L z = Nqe$/\, and in the case of two cylinders 
L z = Nqeo/(2X). We use periodic boundary conditions 
in the z direction by replicating the main simulation box 
infinitely many times in that direction. 

The energy of the system for a given configuration of 
the cylinder (s) is composed of two main parts 



JV 



N 



f3W = J2PW c(r t ) + J2l 3W cc(^,r J ) (Bl) 



i=l 



The first term is the counterion-cylinder interaction en- 
ergy as given by Eqs. ([20)1 and (|34p in the text. The 
second term in (|B1[) is counterion-counterion interaction 
energy which includes also the contribution from interac- 
tions with counterion image charges. Due to the periodic 
boundary conditions used in the simulations, this latter 
term involves infinite summation series (79[. This is be- 
cause the counterions and image charges in the main sim- 
ulation box interact also with their periodic "copies" as 
obtained by the replication of the main simulation box to 
an infinite number of simulation box copies. These sum- 
mations can be evaluated as explained in the forthcoming 
Sections. 

The interaction energy between two given counterions 
i and j can be written as 



oo 



Tj\ < 2R C 



3, 



w (r j , Tj ) + W[ 1U (r j , r j ) otherwise. 



(B2) 

The first equation corresponds to the hard-core repulsion 
between two counterions overlapping counterions. The 
second equation gives the image self-energy of the z-th 
counterion, i.e., it takes into account the interaction of a 
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given counterion with its dielectrically induced images in 
the polarizable dielectric core. The third equation in (|B2[) 
gives the interaction between two different counterions i 
and j. It is written as the sum of two different contri- 
butions which will be calculated below: i) the Coulomb 
interaction, Wo, between the z-th (located at rj) and the 
j-th counterions and all its periodic copies (located at 
Tj + kL z e z for integer k — — oo, 



w {r i ,T j )= 



,00 



k=-c 



r, + kL z e z 



(B3) 



where Aw ff se t accounts for the offset generated at the 
threshold, i.e., 



Awoffsot = w D (Ap = Ap ) - wl{Ap = Ap 



(B8) 



We choose the threshold value as Apj_o = OAL z and 
use a cut off of fco = 6 in the direct scheme and sum up 
to 2 + [1.5/ Ap±] terms in the Lekner scheme, which give 
a relative error smaller than 10 -5 . 



Calculation of Wi n 



and ii) the contribution involving interactions with the 
dielectrically induced image charges Wi m that will be de- 
fined later in this Appendix. 



Calculation of too 



Let us now consider the second contribution that enters 
Eq. (|L32|) . i.e., W[ m (ri,rj), which is the interaction energy 
obtained by summing up all contributions that involve 
dielectrically induced image charges throughout different 
periodic copies of the main simulation box along the z 
axis. It can be written in the case of one cylinder as 



The first contribution wq can be written as 



Ap/L z < Ap , 
Ap/L z > Ap , 



(B4) 



where we have introduced two different summation 
schemes depending on the radial separation Ap between 
the two position vectors and tj in order to increase 
the convergence of the sum over the long-range electro- 
static interactions. If radial distance Ap = \J Ax 1 + Ay 2 
is smaller than a threshold Ap , we first sum up kg terms 
directly and then use an approximation to estimate the 
remaining part of the series (from fco + 1 to infinity). We 
refer to this as the "direct" summation scheme. This 
method is different from the so-called Sperb [HJ sum- 
mation scheme which involves calculating various special 
functions. For our purposes of relative precision up to 
10~ 5 we found it more efficient and can be made in prin- 
ciple as accurate as necessary. It follows as 



w D (rj,rj) 
q 2 £ B 



k— — k( 



n . 2Az 2 -A P 2 ~ 1 



k=k + l 



where the last sum over 1/fc 3 is equal to — ^tp"(ko + 1), 
where ip is digamma function, and needs to be calculated 
only once in the course of the simulations. The term s(k) 
is defined as 



s(k) = 



^jAp 1 + (Az + kL 2 



(B6) 



This scheme gives an error of the order 0{\/k^). 

For larger radial separations (Ap > Apo), we use the 
so-called Lekner summation scheme [8l| as it converges 
more rapidly. It is given by 



2 Ap 4 °° 
= -y- In -J— + y~y~] K (k m Ap) cos{k m Az) 



+Aw ff S ot, 



(B7) 



Wi 



r(Ti, rj) = f3(e q) 2 ^2 u im (r l: r 3 + nL z e z ), (B9) 



n— — oo 



where e z is a unit vector pointing in +z direction. In- 
serting the image kernel (fTTj) into the above equation, 
one ends up with an infinite summation of terms includ- 
ing cos k{Az + nL z ), which can be expressed in terms of 
6 functions, i.e., 



cos k(Az +nL z ) = 7rcos kAz 8{kL z — 2 / Kn). 



n— — oo 



n— — oo 



(BIO) 

This enables us to simply integrate over the wave- vector 
variable k and thus obtain 



2 oo oo 

^ ] ^ ] Cm(k n a)K m (k n pi)K m (k n pj) 

m—O n— — oo 

x cos k n Az cos mAip, (Bll) 



where k n = 2im/ L z . The n = term needs to be handled 
carefully by taking the limit k n — > as 



lim £ m (kna)K m (knPi)K m (knPj) = — (—— 
fc-s-o m \ piPj 



(B12) 



Using the summation identity as used in Eq. (|A9[) . we 
find the final expression 



Wl m (Ti,Tj) 



L, 



1 - 2 



cos Aip 



PiP 3 



+ — ^ ^ tjm,(k n a)K m (k n pi)K m {k n pj) 

rn—O n—1 



x cos k n Az cos mAip. 



(B13) 



This expression is the most computationally expensive 
part of the simulations, and its efficiency decreases for 
Pi/L z <C 1, i.e., for counterions near the cylindrical core. 
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In the simulations we truncate these summations in such 
a way that the relative cut-off error is less than 1CV 4 . 

In the case of two cylinders, we use the first-order im- 
age approximation in the simulations as well, and thus 
the image-charge contributions are simply the sum of the 
contributions for each cylinder separately, i.e., 



,(2) 



i(r<,i' r J,i) + Wim(i\2,r 



3,2), 



(B14) 



where indices 1 and 2 represent coordinates centered at 
the first and the second cylinder, respectively. 



Appendix C: Calculation of the inter-cylinder force 
in the simulations 



In this Section, we explain the method used to directly 
evaluate the force between two charged cylinders within 
MC simulations. The force acting on a given cylinder 
(say cylinder 1) is composed of three contributions, 



TV N 

F = F 00 +J2Foc{i) 

i=l i=l 



^ ^ F Q sm (^) ) 



(CI) 



where i<oo is the bare electrostatic force due to the inter- 
action with the other cylinder, Fo c (i), is the electrostatic 
force due to the interaction of the cylinder with the i-ih 
counterion (summed over all counterions i = 1, . . . , N), 
and -F sm(*) is the force due to the osmotic pressure from 
the z-th counterion. 

The first contribution follows from the interaction en- 
ergy W 00 , Eq. 422), as 



F, 



dW, 



00 



which gives 



00 



00 



2£ 2 



OR 



q 2 £ B R 



1 



2Ac 



R 2 -a? 



(C2) 



(C3) 



The second contribution follows by differentiating Wo c , 
Eq. (|34]l. with respect to R by assuming that the cylinder 



1 and the 7-th counterion are fixed, i.e., 



F 0c (i) 



dW 0c 



dR 



(C4) 



P2,i iP2,i i i P2,i, l P2,i— COnS *' 

which after some algebra gives 



PF Qc (i) 
L 



+a(i 



(1-A) 



R 2 



cos <p itl 

Pi,l 

cos y*i 
Pt,x 



A 



a z cos ipU 



R 2 



Pli 



(C5) 



Here p^\ and tp^i are the radial distance and the az- 
imuthal angle of the i-th counterion position with respect 
to cylinder 1, and p* ix and ip\ x are the corresponding co- 
ordinates with respect to the shifted axis by a 2 /R from 
the axis of cylinder 1 toward the axis of the cylinder 2, 
Fig. [31 Similar definitions apply to p^%, fi t 2, P*2 an d 

* 

The osmotic (third) term in Eq. (jCip results from the 
collisions of counterions with the cylinder. This contri- 
bution amounts to a pressure of counterions exerted on 
the cylinder surface, which is proportional to the contact 
density of counterions at the cylinder surface no. Thus 
the infinitesimal force acting on the cylinder is given by 
dF = jio(ip)kTadip L cos (p. Note that the contact coun- 
terion density no(ip) is given by the number of counteri- 
ons inside a small box of dimensions a Sip x 5p x L next 
to the cylindrical surface. The total osmotic force of all 
ions in simulation cell is then 



= um — y 



cos lfi it i 
~Jp~ 



®{a + 5p-p it r), (C6) 



where 9 is the Hcaviside step function which is 1 only if 
the i-th counterion falls inside a shell thickness 8p from 
the cylinder surface. The above procedure in principle 
gives the exact value of the osmotic force when Sp — > 0. 
It also gives a convenient method to calculate the os- 
motic contribution within MC simulations by evaluating 
the force for different values of the shell thickness and 
estimating the limiting value by extrapolation. 
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